import sympy as sp
import sympy.integrals.quadrature

x,y = sp.symbols('x,y')

expr = sp.exp(-x**2+y)

res1 = sp.integrate(expr,(x,-1,1),(y,-1,1)).evalf()

res2 = 0
n = 4
nodes, weights = sympy.integrals.quadrature.gauss_legendre(n,5)
for i in range(n):
    for j in range(n):
        res2 += weights[i]*weights[j]*expr.subs({x:nodes[i],y:nodes[j]}).evalf()

print(res1,res2)
